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Abstract 

We have developed a general technique to study the dynamics of the quantum adiabatic evolution 
algorithm applied to random combinatorial optimization problems in the asymptotic limit of large 
problem size n. We use as an example the NP-complete Number Partitioning problem and map 
the algorithm dynamics to that of an auxiliary quantum spin glass system with the slowly varying 
Hamiltonian. We use a Green function method to obtain the adiabatic eigenstates and the minimum 
exitation gap, g min = 0 (n 2 ~ n G), corresponding to the exponential complexity of the algorithm 
for Number Partitioning. The key element of the analysis is the conditional energy distribution 
computed for the set of all spin configurations generated from a given (ancestor) configuration by 
simultaneous flipping of a fixed number of spins. For the problem in question this distribution is 
shown to depend on the ancestor spin configuration only via a certain parameter related to the 
energy of the configuration. As the result, the algorithm dynamics can be described in terms of 
one-dimensional quantum diffusion in the energy space. This effect provides a general limitation 
of a quantum adiabatic computation in random optimization problems. Analytical results are in 
agreement with the numerical simulation of the algorithm. 
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I. INTRODUCTION 


Since the discovery by Shor [1] nearly a decade ago of a quantum algorithm for effi- 
cient integer factorization there has been a rapidly growing interest in the development of 
new quantum algorithms capable of solving computational problems that are practically 
intractable on classical computers. Perhaps the most notable example is that of a combina- 
torial optimization problem (COP). In the simplest case the task in COP is to minimize the 
cost function (“energy”) E z defined on a set of 2 n binary strings z = {^ 1 , • • ■ > z n) z j — 
each containing n bits. In quantum computation this cost function corresponds to a Hamil- 
tonian Hp 

^ = X>|z)<z| M 

z 

iz) = ® ^2)2 ® * * ' ® |^n)n- 

where = 0, 1 and the summation is over T states \z) forming the computational basis 
of a quantum computer with n qubits. State \z,)j of the j-th qubit is an eigenstate of the 
Pauli matrix d 2 with eigenvalue Sj = 1 - 2 Zj ( S ) = ±1). It is clear from the above that the 
ground state of H P encodes the solution to the COP with cost function E z . 

COPs have a direct analogy in physics, related to finding ground states of classical spin 
glass models. In the example above bits Zj correspond to Ising spins S r The connection 
between the properties of frustrated disordered systems and the structure of the solution 
space of complex COPs has been noted first by Fu and Anderson [2]. It has been recognized 
[3] that many of the spin glass models are in almost one-to-one correspondence with a 
number of COPs from theoretical computer science that form the so-called NP-complete 
class [4]. This class contains hundreds of the most common computationally hard problems 
encountered in practice, such as constraint satisfaction, traveling salesmen, and integer 
programming. NP-complete problems are characterized in the worst cases by exponential 
scaling of the running time or memory requirements with the problem size n. A special 
property of the class is that any NP-complete problem can be converted into any other NP- 
complete problem in polynomial time on a classical computer; therefore, it is sufficient to 
find a deterministic algorithm that can be guaranteed to solve all instances of just one of the . 
NP-complete problems within a polynomial time bound. It is widely believed, however, that 
such an algorithm does not exist on a classical computer; whether it exists on a quantum 
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computer is one of the central open questions. Ultimately, one can expect that the behavior 
of new quantum algorithms for COPs and their complexity will be closely related to the 

properties of quantum spin glasses. 

Recently, Farhi and co-workers suggested a new quantum algorithm for solving combina- 
torial optimization problems which is based on the properties of quantum adiabatic evolution 
[5], Running of the algorithm for several NP-complete problems has been simulated on a 
classical computer using a large number of randomly generated problem instances that are 
believed to be computationally hard for classical algorithms [6-9]. Results of these numerical 
simulations for relatively small size of the problem instances ( n < 20) suggest a quadratic 
scaling law of the run time of the quantum adiabatic algorithm with n. Furthermore, it was 
shown in [10] that the previous query complexity argument that led to the exponential lower 
bound for unstructured search [11] cannot be used to rule out the polynomial time solution 
of NP-complete Satisfiability problem by the quantum adiabatic algorithm. 

In [10, 12-15] special symmetric cases of COP were considered where symmetry of the 
problem allowed the authors to describe the true asymptotic behavior (n -> oc) of the 
algorithm. In certain examples considered in [5, 13] the quantum adiabatic algorithm finds 
the solution in time polynomial in n while simulated annealing requires exponential time. 
This effect occurs due to the special connectivity properties of the optimization problems that 
lead to the relatively large matrix elements for the spin tunneling in transverse magnetic 
field between different valleys during the quantum adiabatic algorithm. In the examples 
considered in [13] the tunneling matrix element scales polynomially with n. On the other 
hand, in simulated annealing different valleys are connected via classical activation processes 
for spins with probabilities that scale exponentially with n. It was also shown for certain 
simplified examples [14, 15], that quantum adiabatic algorithm can be modified to completely 
suppress the tunneling barriers even if the corresponding classical cost function has local 

minima well separated in the space of spin configurations. 

However, so far there are no study on the true asymptotic behavior of the algorithm 
for the general case of randomly generated hard instances of NP-complete problems. Also 
there are no analysis of the limitations of the quantum adiabatic computation arising from 
the intrinsic properties of disorder and frustration in this problems. Such analysis is of the 
central interest in this paper. 

In Sec. II we introduce the random Number Partitioning problem and describes condi- 
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tional cost distributions (neighborhood properties) in this problem. In Sec. Ill we describe 
the concept of quantum adiabatic computation applied to combinatorial optimization prob- 
lems and introduce a Green function method for the analysis of the minimum gap. In Sec. IV 
we describe the effect of quantum diffusion in the algorithm dynamics, derive the scaling for 
the minimum gap and the complexity of the algorithm for the random Number Partitioning 
problem. We also obtain the scaling of the minimum gap numerically from the form of the 
cumulative density of the adiabatic eigenvalues at the avoided-crossing point. In Sec. V we 
discuss the results of the simulations of the time-dependent Schrodinger equation to sim- 
ulate quantum adiabatic computation for Number Partitioning and obtain its complexity 
numerically. 


II. NUMBER PARTITIONING PROBLEM 


Number Partitioning Problem (NPP) is one of the six basic NP-complete problems that 
are at the heart of the theory of NP-completeness [4]. It can be formulated as a combinatorial 
optimization problem: Given a sequence of positive numbers {ai, . . . , a n } find a partition, 

i.e. two disjoint subsets A and A', such that the residue 


E = 


CL j £.4 


n a J 

ajeA' 


( 2 ) 


is minimized. In NPP we search for the bit strings z = {z u (or corresponding Ising 

spin configurations S = {S u . . . , S„}) that minimize the energy or cost function E z 


n 

E z = |Gs| > Gs = Sj = 1 — 2 Zj, (3) 

i - i 

where Sj = 1 (zj = 0) if a j G A and S 3 = -1 (zj = 1) if % G A'. The partition S 
with minimum residue can also be viewed as the ground state of the Ising spin glass, — 
corresponding to the Mattis-like antiferromagnetic coupling, Jij = —CLidj. 

NPP has many practical applications including multiprocessor scheduling [16], cryptogra- 
phy [17], and others. The best deterministic heuristical algorithm for NPP, the differencing 
method of Karmakar and Karp [18], can find with high probability solutions whose energies 
are of the order l/n aIogn for some a > 0. The interest in NPP also stems from the re- 
markable failure of a standard simulated annealing algorithm for the energy function (3) to 


4 



find good solutions, as compared with the solutions found by deterministic heuristics [19]. 
The apparent reason for this failure is due to the existence of order 2" local minima whose 
energies are of the order of 1/n [20] which undermines the usual strategy of exploring the 

space of the spin configurations S through single spin flips. 

The computational complexity of random instances of NPP depends on the number of bits 
b needed to encode the numbers a,. In what follows we will analyze NPP with independent, 
identically distributed (i.i.d.) random b-bit numbers a r Numerical simulations show [21, 22, 
26] that solution time grows exponentially with n for n « b then decreases steeply for n > b 
(phenomenon of “peaking”) and eventually grows polynomial^ for n > b. The transition 
from the “hard” to computationally “easy” phases at n * b has features somewhat similar 
to phase transitions in physical systems [23]. The detailed theory of the phase transition in 
NPP was given in Refs. [24, 25], If one keeps the parameter £ = b/n fixed and lets n -4 oo 
then instances of NPP corresponding to <; > 1 will have no perfect partitions with high 
probability. On the other hand for f < 1 number of perfect partitions will grow exponentially 
with n. Transitions of this kind were observed in various NP-complete problems [28]. In 
what follows we will focus on the computationally hard regime f > 1- 

A. Distribution of signed partition residues 

In Fig. 1 we plot an array of 2 n partition energies E z = |O z | sorted in increasing order. 
While the values of individual energies are random and depend on the particular instance 
of NPP (i.e., the set of numbers a 3 ) it can be inferred from Fig. 1 that on a coarse- 
grained scale (i.e. after averaging over individual energy separations) the form of the typical 
energy distribution is described by some universal function for randomly generated problem 
instances. To describe it we introduce for a given set of randomly sampled numbers a j a 
coarse-grained distribution function of signed partition residues (3) 

P(ft) = 2- n — / dr i 6 ^ 

AS2 Jq-AQ/ 2 z€{0,l} n 

Here <S(x) is the Dirac delta-function; the sum is over 2" bit-strings z and 2'” is a normal- 
ization factor. In (4) we average over an interval AS2 of the partition residues whose size is 
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( 5 ) 


chosen self-consistently, AQ 2 n /P(Q). Using (3) we can rewrite (4) in the form 


p(n) = hT ds< (^r) /w cos(ns) ’ 

n 

1(8) = PJcos(Gjs), Q(x) = sin(x)/x. 
j = 1 

Here ((x) is a window function that imposes a cut-off in the integral (5) at s ~ 2/Afi. For 
large n this integral can be evaluated using the steepest descent method. In the following 
we shall assume that the 6-bit numbers aj are distributed inside of the unit interval [0, 1] 
and are integer multiples of 2~ b , the smallest number that can be represented with available 
number of bits 6. We note that for large n the function I(s) has sharp maxima (minima) 
with width ~ n -1 ^ 2 at the points s* = ki\2 b , k — 0,1....; |/(s*;)| = 1. Only one saddle 
point at s = 0 contributes to the integral in (5) due to coarse-graining of the distribution (4). 
Indeed, it will be seen below that the window size 2/Af2 can be chosen to obey the conditions 
1 < n 1/2 /Afl <C 2". Therefore in the case of high-precision numbers, 6 » n, saddle-points 
Sjt with k > 0 lie far outside the window and their contributions can be neglected (see also 
Appendix A). On the other hand the window function ((x) can be replaced by unity while 
computing the contribution from the saddle-point at s = 0. Finally we obtain for |fi| n 
(cf. [29]) 


P(Q) = 


a 2 (0) = 


y/2 7t<7 2 (0) n 


exp 


Q 2 


2cr 2 (0) n 


+ 0(n- 3/2 ) 


n 




j = i 


( 6 ) 


The coarse-grained distribution P(£l) depends on the set of a^’s through a single self- 
averaging quantity <r(0) (cf. [23]). 

One can also introduce the distribution P(E) of cost values (energies) E x = |fi,[. Due to 
the obvious symmetry of the NPP, the cost function E z in (3) does not change after flipping 
signs of all spins, Sj — > —Sj. Therefore 


P(E) = 1/2 P(±E). 


( 7 ) 


We emphasize that, according to Eq. (6) for a typical set of high-precision numbers a 3 the 
energy spectrum in NPP is quasi-continuous, and there are only two scales present in the 
distribution P(E): one is a “microscopic” scale given by the characteristic separation of the 
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• • • 7Ti cinrl snnthpr is ffiven bv the mean partition energv (E) 

individual partition energies, E m in , and another is given uv 

(or the distribution width ( E 2 ) 1/2 ) 

r(Q) n l < 2 2- n , (E 2 ) = |(£) 2 = na 2 ( 0). (8) 


^min ^ 


This justifies the choice for Afl above that corresponds to coarse-graining over many indr 
vidual energy level separations. 



FIG. 1: Dotted plot is a normalized distribution of partition energies E, = |fi.l defined in 
Eqs. (4), (7). Size of the input n = 20, precision 6 = 35. Insert: plot of partition energies E , 
sorted in increasing order (fc gives the position of partition energies in a sorted array). Asymp- 
totic result based on Eq. (6) is visually indistinguishable from the exact result. 

We note that the distribution P{ fl) (6) is Gaussian for E A' n and can be unde 
in terms of a random walk with coordinate SI using Eq. (3). The walk begins at the origin, 
n = 0, and makes a total of n steps. At the j-th step SI moves to the right or to the left 
by “distance” 2 o, if S, = 1 or S, = -1, respectively. In the asymptotic limit of large n the 
result (6) corresponds to equal probabilities of right and left moves and the distribution of 

step lengths coinciding with that of the set of numbers {2 a^}. 

Finally, the energy distribution function P(E) of the form (6), (7) was previously obtained 
by Mertens [29] using explicit averaging over the random instances of NPP. He also computed 
the partition function Z(T) for a given instance of NPP at a small finite temperature T using 
the steepest-descent method and summation over the saddle-points s t = kx 2 1 similar to our 
discussion above [23] (in his analysis k B T played a role similar to our regularization factor 

AQ in (4), (5)). 

We emphasize however, that the approach in Ref. [23] based on Z(T) is necessan y 
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restricted to the analysis of the “static” properties of NPP at E ~ 2“ n , i.e., the phase 
transition in the number of perfect partitions [23] when the control parameter f = njb 
crosses a critical value. On the other hand distribution P($l) (4) introduces at finite energies, 
as well as the conditional distribution introduced in the next section also allow us to directly 
study the intrinsic dynamical properties of the problem in question such as the dynamics of 
its quantum optimization algorithms. 


B. Conditional distribution of signed partition residues 


Consider the set of bit-strings z' obtained from a given string z by flipping r bits. The 
conditional distribution of the partition residues Vl z > (3) in the r-neighborhood of z can be 
characterized by its moments: 

<Q fc > = (”) 1 E ^ W,z)> k = 1,2,... (9) 

' ' z , e{o,i}’ 1 

Here <S m>t is a Kronecker delta and function D{ z,z') computes the number of bits that take 
different values in the bit-strings z and z'. It is the so-called Hamming distance between the 


two strings 

D(z,z’) = EN ~ z i\- ^ 

The Hamming distance v = D( z, z ; ) between the bit-strings is directly related to the o\erlap 
factor q between the corresponding spin configurations often used in the theory of spin 


glasses [3, 29]: 


= ~YS J S’ = 1 - -D{ z,z'). 
n ^ J n 


(ii) 


(in what follows we shall use both quantities r and q ). For k = 1, 2 in (9) one obtains after 
straightforward calculation the first and second moments of the conditional distribution 


(0) = q£l z, 


1 fiz \ 

n (E 2 )J ’ 

(12) 

(ft 2 ) - (f l } 2 = na 2 (q) | 

! + n-l) 0 

(13) 

a (q) = cr(0) (1 - q 2 ) 1/2 , 

2r 

q = 1 , 

n 


(14) 


where <r(0) and (E 2 ) are given in (6) and (8), respectively. 
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The conditional distribution of ft z » can also be defined in a way similar to (4) 


/ \ —1 -i + /2 ^ 

P rz (0')=f n ) — - / dr] 22 5{v - <5r,o(z',z) 

’ W in'-AO'/2 z'e{o,i} n 


( 15 ) 


z'e{o,i}" 

where averaging is over the small interval AST that, however, includes many individual 
values of n, for a given r. It is clear from (12), (13) that the first two moments of P v (tt) 
depend on z only via the value of Q.. This does not hold true, however, for the higher- 
order moments that depend on other functions of z as well. For example, {Q 3 ) involves the 

quantity e ^ c * 

Our main observation is that in the asymptotic limit of large n the conditional distri- 
bution P rz (Q') is well-described by the first two moments (12), (13). Then, according to 
the discussion above, its dependence on z is only via ft.. The detailed study of the higher 
moments (9) will be done elsewhere. Here we use the following intuitive approach relevant 
for analysis of the computational complexity of the quantum adiabatic algorithm for the 
NPP . We average P r , z (fl') over the strings z with residues ft z inside a small interval All 
(containing, however, many levels ft.). After such averaging the result, P r {tt |ft), can be 
written in the form 


Pr(ft'lft) = 


PrjO^ 

P{n) ’ 




osrt+An/2 


dr] E %-ft z )P r)Z (ft'), 

/n- An/2 z6 {o,i} n 


(16) 

(17) 


where P(Q) is given in (6). We note that Eq. (16) formally coincides with the Bayesian 
rule expressing the conditional distribution function P r (tt\Q) through the 2-point (joint) 
distribution function P r (ft',ft) and the single-point distribution P(Q) (6). 

The explicit form of P r (fl|ft) is derived in Appendix B in a manner similar to the deriva- 
tion of P(E) in Sec. II A. The results are presented in Eqs. (B9) and (BIO). They show 
that P r (ft', ft) in the limit n » 1 is indeed well described by its first two moments that cor- 
responds precisely to the expressions given in Eqs. (12), (13) above. From this we conclude 


that 

p r , z (tf) = PrW in.). (18) 

In the case r = 1 there are n strings z' at a Hamming distance 1 from the string z. 
Partition energies corresponding to these strings equal |fi s — 2 o j^I> 1 — 2 — n ( c ^' ^)) 
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After the coarse-graining over the energy scale 0(l/n) in the range, |fi|, Iff I « «. *>*' 
conditional distribution z is a step function in the interval 0, - ST 6 [-2, 2], For r = n- 1 
one has the same form of the distribution but for + fi'. Both results correspond to 
nearly equal distribution of spins between between ±1 values. Then in the range of energies 

< 1 o ne has: 

P r<z (rt) ~ P r = 1/2 + O ^ , r = 1, n — 1 (n>l) (19) 

For r,n - r > 1 distribution P r , z (ft') has a Gaussian form with a broad maximum at 
O' = q n z (cf. Eqs. (12),(13),(B9)) . Near the maximum we have: 


P rz { O') « P T = . l l°'l’ 1^1 ^ ri 1/2 a{q). (20) 

5 ^J2'Kno 1 {q) 

We studied the conditional distribution in NPP numerically as well (see Fig.2 and Sec.B). 
The results are in good agreement with theory even for modest values of n < 30. 

The characteristic spacing between the values of the partition residues in the subset of 


strings t! with D( z', z) = r is l/(P r (”)) for not too large E z , E z , (see above). This spacing 
decreases exponentially with the magnitude of the string overlap factor, \q\ = \(n- 2r)/n\. 
The hierarchy of the subsets corresponding to different values of \q\ form a specific structure 
of NPP. We note that the distribution of partition residues within the hierarchy is nearly 


independent of the ancestor string z in a broad range of energies E’ < n l < 2 where P r ,z{E') ~ 
P T . One can see that the magnitude of the overlap factor q between two strings with energies 
within a given interval [0 ,E) is limited by some typical value q satisfying the following 


equation: , _ 

®(")ft = 1’ W = l ~ 2V n <21) 

The smaller E is, the smaller |g| is: strings that are close in energy are far away in the 
configuration space. This property gives rise to an exponentially large number of local 
minima for small values of E , that are far apart in the configuration space. For example, 
strings with E z ~ B mln typically correspond to |,i = 0( 1/n), they can be obtained from 

each other only by simultaneously flipping clusters with - n/2 spins. 

Eq. (21) describes the dynamics of a local search heuristic (e.g., simulated annealing). 
It shows that the average cost value E during the search decreases no faster than 0(1/M) 
where M = O ((")) is the number of generated configurations. This result coincides with 
that obtained in [29] using a different approach. It says that any classical local search 
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FIG. 2: Plots of the (scaled) conditional distribution (15) a = cr(0)(27m) 1/2 (AQ) 1 f 0 dr] P r , z {r]) 
vs r are shown with points. We use coarse-graining window AQ=0.3. Different plots correspond to 
29 randomly selected bit-strings z with energies |fi,| £ [0, 0.3] for one randomly generated instance 
of NPP with n = 30 and b = 35. For r, n - r > 1 the values of s corresponding to different strings 
axe visually indistinguishable from each other. Dashed line is a plot of a(0)/o{q) vs r given in 
(13) {q = 1 - 2 r/n). Insert: plots of the integrated quantity given in (Bll), Q = \ J 0 dr] P r<z (v) 
vs x = n/(a(qW 2^), for different values of r = 2, . . . , n/2 and randomly selected bit-string z with 
energy |fi,| close to 0. All plots correspond to the same instance of NPP as the main figure. Plots 
for different values of r are visually indistinguishable from each other and from the theoretical 

curve given in (B12). 

heuristic for NPP cannot be faster than random search. Indeed, during local search the 
information about the “current” string z with E z < 1 is being lost, on average, after one 
spin flip (cf. Eqs. (19), (20)). We show below that precisely this property of NPP also leads 
to the complexity of the quantum adiabatic algorithm corresponding to that of a quantum 

random search. 

We note that one can trivially break the symmetry of NPP mentioned above by introduc- 
ing an extra number a„ and placing it, say, in the subset A. In this case different partition 
energies will still be encoded by spin configurations S = {Sr, . . . , S„} (or corresponding bit- 
strings z) with fl s = oo + Sr a i “d E. = |f) s | (cf. 3). We shall adopt this approach in 
the analysis of the performance of the quantum adiabatic algorithm for NPP given below. 
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III. QUANTUM ADIABATIC EVOLUTION ALGORITHM 


In the quantum adiabatic algorithm [5] one specifies the time-dependent Hamiltonian 
H(t) = H(t/T ) 

H(t) = (1 -r) V + tHp, ( 22 ) 

where T = t /T is dimensionless “time”. This Hamiltonian guides the quantum evolution 
of the state vector \^(t)) according to the Schrodinger equation id\ip{t))dt = H(t)\i>(t)) 
from t = 0 to t = T, the run time of the algorithm (we let h = 1). H P is the “problem” 
Hamiltonian given in (1). V is a “driver” Hamiltonian, that is designed to cause transitions 
between the eigenstates of H P . In this algorithm one prepares the initial state of the system 
ip(0) to be the ground state of H(0) = V. In the simplest case 

v = -iy„ wo)) = 2-/’D*>. (23) 

j = i z 

where is a Pauli matrix for j-th qubit. Consider instantaneous eigenstates | M T )) of 
H(t) with energies A^(r) arranged in nondecreasing order at any value of r G (0, 1) 

H\4>r,) = K\<l>n), V = 0,l,...,2 n -l. (24) 

Provided the value of T is large enough and there is a finite gap for all t € (0 ,T) between 
the ground and excited state energies, g(r) = Ai(r) - A 0 (r) > 0, quantum evolution is 
adiabatic and the state of the system | ip{t)) stays close to an instantaneous ground state, 

| <f> 0 (t/T)) (up to a phase factor). Because H(T) = H P the final state \ip(T)) is close to the 
ground state |0 o (r = 1)) of the problem Hamiltonian. Therefore a measurement performed 
on the quantum computer at t = T (r = 1) will find one of the solutions of COP with large 
probability. 

There is a broad class of COPs from theoretical Computer Science where the number of 
distinct values of a cost function scales polynomially in the size of an input n. An example 
is the Satisfiability problem in which the cost E z of a given string z equals the number 
of constrains violated by the string. For those problems, the spectrum of H(t), at the 
beginning (r « 0) and at the end (r » 1) of the algorithm, consists of a polynomial number 
of well-separated energy levels. Quantum transitions away from the adiabatic ground state 
occur most likely near the avoided-crossing points r « t* where the energy gap g{r) reaches 
its minima [9]. Near the avoided-crossing points, the spectrum of H(t) is quasi-continuous, 
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with the separation between individuals eigenvalues scaled down with n. The probability of 
a quantum transition, 1 - \m)\Mt/T))\l T , small provided that 

T>> m±L Mk=l, 9m , n = o min [A.M-AoW], (25) 

9 min 

{H t = dH/dr). The fraction in (25) gives an estimate for the required runtime of the 
algorithm and the task is to find its asymptotic behavior in the limit of large n » 1. The 
numerator in (25) is less than the largest eigenvalue of H T = H P - V , typically polynomial 
in n [5], However, g min can scale down exponentially with n and in such cases the runtime 
of the quantum adiabatic algorithm will grow exponentially with the size of COP. 

A. Implementation of QAA for NPP 

As suggested in [5] the quantum adiabatic algorithm can be recast within the conventional 
quantum computing paradigm using the technique introduced by Lloyd [30). Continuous- 
time quantum evolution can be approximated by a time-ordered product of unitary op- 
erators e - i ( 1 ~~ r i‘) vs e ~iTkHp6^ corresponding to small time intervals {tk,tk + ^)- Operator 
e -i(i-T k )V6 typically corresponds to a sequence of 1- or 2-qubit gates (cf. (23)). Operator 
e -ir k H P 6 is diag onal in the computational basis |z) and corresponds to phase rotations by 
angles EJ. Since in the case n « b , the average separation between the neighboring values 
of E x is 1 jP[E) = 0{ 2~ n ), the quantum device would need to support a very high precision 
in its physical parameters (like external fields, etc.) to control small 0(2-) differences in 
phases. Since this precision scales with n exponentially it would strongly restrict the size 
of an instance of NPP that could be solved on such a quantum computer. This technical 
restriction is generic for COPs that involve a quasi-continuous spectrum of cost-function 
values. Among the other examples are many Ising spin glass models in physics (e.g., the 
Sherrington-Kirkpatrick model [3]). To avoid this restriction we introduce a new oracle-type 
cost function £ z that returns a set of values 

£ z = c(f2 z ), c{x) — > {eo,Si,...,e M } (e fc+ i > e k ), ( 26 ) 

that can be stored using a relatively small number of bits <9(logn). For example, we can 
divide an interval of partition energies (0 ,B), B = £J=o% into bins whose sizes grow 
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exponentially with the energy. Then the new cost will take one value per bin 


c{ X ) = £k = ~M + k for u k < W < uJk+i, 

u k = (2 k - 1) A, k = 0, . . . , M. 


( 27 ; 


The last bin is co M <\£l z \<B where we have £ z = £m - 0- The value of the cutoff u M < B 
is discussed below. In this example the Hilbert space of 2" states |z> is divided into M + 1 
subspaces £ fc , each determined by Eq. (27) for a given k 

Hp = E £ ‘EI Z )< Z I- (28> 

A:=0 zeCfc 

Note that subspace Co contaius the solution(s) to NPP. Dimension d 0 of Co is controlled 
by the value of A in (27) which is another control parameter of the algorithm. We set 
A = 2~ n K/P(0) where the integer /V zs ri. > i is independent of n and determines how 
many times on average one needs to repeat the quantum algorithm in order to obtain the 
solution to NPP with probability close to 1. 

Operator H e projects any state | </>} onto the states with partition residues in the range 
0 < |fi z | < If we choose 

1 <««<(£), (29) 

then the distribution function (6) is nearly uniform for |Q Z | < u M . Therefore the dimensions 
of the subspaces £ fc grow exponentially with t. d k = do 2 fc for fc < M. This simplification 
would not affect the complexity of a quantum algorithm that spends most of its time m ‘-an- 
nealing” the system to much smaller partition residues, ui.m » \Q,\ ~ £min - 0{n 2 ). 

We note that the new discrete- valued cost function defined in (27) is non local. Unlike 
problems such as Satisfiability, it cannot be represented by a sum of terms each involving a 
small number of bits. To implement a unitary operator with H P given in (28) one 

needs to implement the following classical function on a quantum computer 


£ z = 0(ojm — |n.|) 


f A 4- 

l0& l aT^ 


= 2z } ). 

j = i 


(30) 


Here [x] denotes the integer part of a number x\ 0(x) is the theta-function (0(x) 1 for 

x > 0 and 0(x) = 0 for x < 0). The implementation of (30) with quantum circuits involves, 
among other things, the addition of n numbers together with their signs to compute Q z , and 
taking the discrete logarithm of a 6-bit number with respect to base 2. These operations 
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can be performed using a number of quantum gates that is only polynomial in n and b (cf. 

[1] for the implementation of the discrete logarithm). 

Since the implementation of a cost function (26), (30) does not add an exponential over- 
head to the complexity of QAA the feasibility of this algorithm for NPP depends on the 
scaling of the minimum gap g m \ n with n. 


B. Stationary Schrodinger equation for adiabatic eigenstates 


We now solve the stationary Schrodinger equation (24) and obtain the minimum gap g m in 
(25) in the asymptotic limit n -> oo. To proceed we need to introduce a new basis of states 
|x) = \xi)i ® 1 3 : 2)2 ® • • ' ® kn)n where state | xj)j is an eigenstate of the Pauli matrix o x for 
the j-th qubit with eigenvalue 1 - 2x } = ±1. Driver Hamiltonian V can be written in the 

following form: 

v = £v„r", r"= E MM- (31) 

m= 0 * 1 4 h x n =m 

For a particular case given in Eq. (23) we have V m = 2 m-n. Matrix elements of I™ in a 
basis of states |z) depend only on the Hamming distance D( z, z'j between the strings z and 


Wry) = /?,„■), 


c = 2 _n E £ 

q—0 p = 0 

We now rewrite Eq. (24) in the form 


n — r\ r 



( 1) P A m ,q+p- 


(32) 

(33) 


|<A) = . T — u H P W), a = a(r) = 1 - t, (34) 

A — atv 

(we drop the subscript r) indicating the number of a quantum state and also the argument r 
in 4> and A). From (27)-(34) we obtain the equation for the amplitudes 4> z = (z| <f>) in terms 
of the coefficients 

[1 — tG 0 c(fi z )] (f> z = t 77 + T ^ Gd(z,z') <t> z' c(D z ), (35) 

A “ aV ° z'*z 


$ = c(p z ') $ Z ', 


jm 
1 r 


G ' s = E 0£r£ ”' 


m=l 
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Here we separated out a “symmetric” term oc 2 n 4> corresponding to the coupling between 
the states |z) via the projection operator 1° (31). 


IV. MINIMUM GAP ANALYSIS 


A. Coarse-graining of the transition matrix 


We now make a key observation that <j> z in (35) can be determined based on the properties 
of the conditional distribution P T , X (E) (15) and the form of the Green function G r ( A). We 
sum the Green function G d(m -) over all possible transitions from a given state z' to states 
z' ^ z with energy e k ■ For not too large partition residues of the initial and final states we 

obtain 


Gd{z : z ')(^0 ~ fz',k{^) 

zeC k ,z^z r 


Ft{ A) 


MfC A) 

2 m- 




dr 


a(0) /n 


a (1 — 2r/n) \r 


G r (A) 


(36) 

(37) 


|fi B | < (£), 


2u)m 

^AEY 


(38) 


Function <j(g) above is defined in (14) and / z ',*( A) is a small correction described below. In 
function s(A) we replaced summation over the integer values of r by an integral. It can be 
evaluated using the explicit form of G r (A) that decays rapidly with r. In what follows we 
will be interested in the region |A — Q:Vo| < 1 where 

—2a G r (A) =(")"’ E _ 2 -” (In r + 7 ) . (39) 


(7 is Euler’s constant) and s(A) ~ — ln2/(2o:). We note that 


-2aG r (X) « - 


(n/2 — r) 


, n/2 - r » 1. 


(40) 


Therefore the integrand in s( A) is a smooth function of r for r < n/2 and quickly decays to 
zero for r > n/2. The contribution to the integral in s(A) from the range of r <C n is small 

(0((r/n) 1/2 ). 

We note that term F k in (36) provides an “entropic” contribution to the sum in (36). It 
comes from the large number of states z £ C k corresponding to large Hamming distances 
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r from the state z\ 1 « r < n/2. Each state contributes a small weight, G r oc (?) , and 

number of states for a given r is large, (w»+i - at) © ft > 1- Here ("*+■ ' "* 1 1S an 
energy bin for the subspace C t and ft is the conditional density of states described in Sec. 
II The size of the bin scales down exponentially with k (cf. (27)) and so does the entropic 
term Ft- Below a certain cross-over value of k one has |F*| <£ |/z',fc(A)l- In this case the 
dominant contribution to the sum (36) comes from the states z with small r = D( z, z') ~ 1. 
In particular for k = 0 one can obtain 


/*», o(A) ~ Gi(X) ^2 + G>( 


where the higher-order terms correspond to D( z',w) > 2. According to (39), |G,(A)| ~ n" 2 
and therefore |/,.„| is exponentially larger than the entropic term, |F„| ~ ~ do 2"". We 

note that, unlike the entropic term, /,,o strongly depends on z' due to the discreteness of 
the partition energy spectrum (<v„n < 1). E.g., depending on a state z', in this case there 
could be either one or none of the states w £ £„ in the sum (41) satisfying D( z'.w) = 1. 


B. Extended and localized eigenstates 

Based on the discussion above we look for solution of Eq. (35) m the following form. 

<j> z = + u z , Z i C Q , ( 42 ) 

where we have explicitly separated out a part of the wavefunction that depends on z 

only via the corresponding value of the partition residue. It satisfies the following equations: 


[1 — t G o(A) c(O)] f(f2) — 


A — a Vo 


j 

+ r / dav(n')c(n')x(n',n,\), 


x(^'> n, a) = ^2 (”) Gr ^ p r( Q '\ty- 

r= 1 ' ' 


where $ is given in (35) and function c(x) takes a set of discrete values (26). Using (35), (42) 
and (43) we obtain equations for u z 


m 

fl -TG 0 {\)e k ] u z = T J2 £ k’ ^2 G D(z,z')(A)u z - +T £ 0 2^ G D{ z,w)(A)</>v 


z e Ck- 


k ’= 1 z 
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Decomposition (42) is only applied to amplitudes <p z with z ^ Cq. The system of equations 
for the components u(Q) and u z is closed by adding Eq. (35) for the amplitudes </> w with 
w G £ 0 (ground states of the final Hamiltonian Hp ) and taking (42) into account. We note 
that Eq.(43) for v(fl) is coupled to the rest of the equations only via the symmetric term <f> 

$ = $ + $ + $ 0 (46) 

/ OO 

dx P(x) v(x) c(x), (47) 

-OO 

M 

$ = ^ ^ ^ ^ U z , $0 = £() ^ 

fe=i z eCk wec 0 

where distribution P(Q) is given in (6). 


1. Minimum gap estimate for ljm ^ (P) 


We will analyze the above system of equations (42)-(47) assuming that the cutoff fre- 
quency u) M satisfies Eq.(29). This condition corresponds to the linear region in the plot of 
the cumulative density of states given in insert to the Fig. II A. According to Eqs. (6), (19), 
in this range the distribution functions P(f2) ~ P„/ 2 and P r (0 , |D) ~ P T take nearly constant 
values and spectral function x(f^> equals 


X (D',D,A) 


s( A) 

yj2ix ncr 2 (0) 


(48) 


where s(A) is given in (37). In this approximation, we can compute $ using equations for 
u z in (45) and also the relations in (36), (37) 


$ = -k(t[is(\)) < f > 0 , K ( x ^ = T+x' ^ 

In the initial stage of the algorithm the amplitudes <f) w of the “solution” states are small 
|$ 0 | = 0{ 2~ n / 2 ). According to (49), we also have |$| = 0{2~ n l 2 ). Neglecting these terms 
and setting <$ « <5, Eq. (43) gives a closed-form algebraic equation for A 

1 + 2 rp, ( — — + •s(A)^ = 0. (50) 

Expanding in a small parameter p, < 1 (cf.(29),(38)), we obtain the eigenvalue 

Aq(t) o(r) V 0 - 2 r/i - — ^ \- 0(p?) (a » /x), (51) 
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that accurately tracks the adiabatic ground state energy, A 0 (r), from r = 0, up until small 

vicinity of the avoided-crossing, r « r* (see below) where |$o| ~ 1- 

In the avoided-crossing region, branch A ‘(r) intersects with another branch, X( [(r), that 
tracks A 0 (r) in the interval of time following the avoided-crossing, r < r < 1. This branch 
corresponds to $ < <f*o, $• It can be obtained from simultaneous solution of equations for u z 
(45) and K that are approximately decoupled from Eq. (43) after * is neglected. Keeping 
this term in (45) gives rise to repulsion between branches A ^(r) at r = r* that determines 

the minimum gap <? m in ( see below). 

To proceed, we obtain the equation for 4> 0 by adding equations for amplitudes K that 
correspond to different states w € £ 0 and neglecting the coupling between these states 
separated by large Hamming distances, D{ w,w') ~ n/2. It can be shown using Eqs. (35) 
and (41)-(45) that u z enters equation for 4>o through the term 


t 2 £ o ^2 o(A)«z, ' ' 

z&Co 

which is is a self-energy term corresponding to elementary bit-flip processes with initial and 

final states belonging to the subspace Co (loop diagrams). 

To express a, in (52) through we solve Eq. (45) using order-by-order expansion in a 
small parameter n~ l (cf. Eqs. (36)-(41) and discussion there). In particular, one can show 
that to the leading order in n~' the self-energy term (52) is determined by lowest-order loops 
with two bit flips that begin and end at £„. Then after some transformations, the equation 


for $ 0 takes the form 


$0 A - T€q - 


TCr£ 0 


E 




A jri X-t£ z 

z'^Cq 


= \£ 0 Td 0 2~ 


4> 


A — aVo 


4- 4> s(A) 


(53) 


Here a = 1 - t (cf. (34) and % is defined above. We now solve Eq. (53) jointly with (43) and 
obtain a closed-form equation for A. We give it below in the region of interest |x - 1/2| « 1 


(A - AIM) (A - AiM) = -n 2 2-"A74 

A « 4 /2 (1 + [IT* In 2 + 0(fi 2 )) , 


where the branch Aj,(r) is given above and the branch A f 0 (r) satisfies Eq. 
there set to zero, 

A o(t) « T£ 0 - 1/2, | t - 1/2| < 1- 


(53) with r.h.s. 

(55) 


19 



(56) 


Avoided-crossing in (54) takes place at r = r* 

Ai< T *) = A'(r-), r-sA + i-log^. 

The value of minimum gap between the two roots of (54) equals 

gmm = nA 2 _n / 2 . (57) 

where A is defined in (27). 

Based on the above analysis one can also estimate the matrix element |{<Pi|# t |<Ao)|t=t* ~ 
n. Then from Eq. (25) (see also discussion after Eq. (28)) one can estimate the run-time of 
the quantum adiabatic algorithm 

T » ^JA»li = 0((ndo)- 1 2"). (58) 

9min 

It follows from the above that eigenvalue branch A{>(t) corresponds to a state, 

ze{o,i} n 

which is extended in the space of the bit configurations |z): according to (43) it contains a 
large number (0(2”)) of exponentially small (0( 2~ n / 2 )) individual amplitudes. This state 
originates at r = 0 from the totally symmetric initial state |tHO)) (23). In the small region 
\ T - r*| ~ g min it is transformed into the state that corresponds to the eigenvalue branch 
Ao(r) and is localized in Hamming distances D( z, w) near the subspace w G Co containing 
the solution to NPP |0 O ) ~ Ew 6 £ 0 l w )- Minimum S a P at the avoided-crossing is determined 
by the overlap between the extended and localized states. 

At later times r > r* a similar picture applies to the avoided crossing of the extended- 
state energy Aj,(r) with energies of localized states A£(r) corresponding to z G £* with 
1 < k < n (excited levels of the final Hamiltonian H P (28)). The existence of the extended 
eigenstate of H(t) whose properties do not depend on a particular instance of NPP follows 
directly from Eq. (43) that involves only a self- averaging quantity x(^> ^ ^)- This quantity 
varies smoothly over the broad range of partition residues |fT|, |fi| < (E) and does not 
allow for the compression of the wave-packet v(fi*) on the much smaller scale 0( 2 _n ). This 
gives rise to an eigenstate with probability amplitude of individual states |z) that depends 
smoothly on energy in this range. 
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FIG. 3: Dots correspond to the plot of the ground state amplitude (z#o) vs partition residue 
|fi z | evaluated at the avoided crossing point r = t* (thin lines connecting the dots are for display 
purposes). Simulations are done for the randomly sampled instance of NPP with n - 10 and 
b = 20; the corresponding value of r* « 0.5. In simulations we relax the condition (29) and the 
value of M in (27) is set automatically to be an integer closest to log 2 Yhj=o a i ( c ^' (^7))- 
Dotted curves are the plots of the two lowest eigenvalues of H(r) vs r for the same instance of 
NPP as in the main figure. Solid lines that start at r = 0 correspond to A = (1 - r)n + fc with 
k _ o, 1 (cf. (51)). Solid lines that ends at r = 1 correspond to A = r e k with k = 0, 1 (cf. (55)). 

2. Analysis of the general case 

The above picture of avoided-crossing remains qualitatively the same when the condition 
(29) is relaxed (cf. insert in the Fig. 3). Away from the avoided-crossing point, r < r*, 
the ground state wavefunction v(Q z ) and energy A * (r) are obtained directly from Eq. (43) 
with replacement * * and Eq. (44) taken into account. Because the spectral function 
x(fi, Q', A) changes only slightly on the scale E m m = 0(n ^ 2 ) the wave packet v ( z ) I ) 

remains extended, |u(fi z | = 0(2 n / 2 ), and therefore $o — 0(2 ). 

Beyond the avoided-crossing point, r > r*, the ground state is localized near w and 
eigenvalue branch A {(r) is obtained from Eq. (53) with r.h.s. set to zero (cf. Sec. IV Bl). 
The point r = r* is located at the intersection of the two branches A£,(r) » X f 0 (r) and the 
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level repulsion is of the order of the overlap factor between the extended and localized states 

5mln ~ £ „(£).) ~ (59) 

w€jCo 

Ground-state wavefunction cp z at the avoided-crossing is shown in Fig. 3 for modest value 
of n, but the separation into slowly- and rapidly-varying parts (42) is clearly seen. 

We did not perform a direct numerical study of the dependence of g m [ n on n since we 
only simulated adiabatic eigenvalues for small instances of NPP. We argue, however, that 
even for a fixed n the scaling of g m in with n can be inferred from the shape of the cumulative 
density of states 

r A 

T}{\)= dxY, 6 (\ k -x), k m = 2"-l, (60) 

k=0 

where X k = X k {r) are eigenvalues of H{r) (24). These eigenvalues are plotted in Fig. 4 near 
the avoided-crossing t = t* where the spectrum of A*, is quasi-continuous. The shape of the 
plot is well approximated by the square-root function: 

A^ « const + > Vm = 2 n ). (61) 

It is clear that for g » 1 we have A,, ~ 2“ n/2 which corresponds to Eq. (57). Note that 
this qualitative analysis is based on the assumption that the asymptotic properties of A 0 for 
large n can be inferred from the behavior of X v for ?? A§> 1 - 

V. SIMULATIONS OF TIME-DEPENDENT SCHRODINGER EQUATION 

We also study the complexity of the algorithm by numerical integration of the time- 
dependent Schrodinger equation with Hamiltonian H(t) and initial state |^(0)} defined 
in Eqs. (22), (23), (27), (28). Here we relax the condition uj m < (E) used above in the 
analytical treatment of the problem; in simulations the value of M is set automatically to 
be an integer closest to log 2 £” =0 ^ (cf. (27)). We introduce a complexity metric for the 
algorithm, C(T) = (1 + T)d 0 /p 0 (T) where p 0 (t) = Ew e c 0 llM*)| 2 . A typical plot of C(T) 
for an instance of the problem with n - 15 numbers is shown in the insert of Fig. 4. At 
very small T the wavefunction is close to the symmetric initial state and the complexity is 
~ 2 n . The extremely sharp decrease in C(T ) with T is due to the buildup of the population 
Po(T) in the ground level, £ z = £o> as quantum evolution approaches the adiabatic limit. At 
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FIG. 4: Dotted line is plot of \ v vs rj at the avoided-crossing r = t* . It is obtained from the 
numerical solution of the stationary Schrodinger equation for the same instance of NPP as in 
Fig. 3. Solid line is a square-root fit A = —6.3 + 0.35 77 1 / 2 (solid line is almost undistinguishable 
from the dotted line). 

certain T = T min the function C(T) goes through the minimum: for T > T min the decrease 
in the number of trials cIq/po(T ) does not compensate anymore for the overall increase in 
the runtime T for each trial. For a given problem instance the “minimum” complexity 
C'min = C(T min ) is obtained via one dimensional minimization over T. The plot of the 
complexity C m \ n for different values of n in Fig. 1 appears to indicate the exponential 
scaling law, <7 min ~ 2° 8n for not too small values of n> 11. 

VI. DISCUSSION 

In conclusion, we have developed a general method for the analysis of avoided-crossing 
phenomenon in quantum spin-glass problems and used it to study the performance of the 
quantum adiabatic evolution algorithm on random instances of the Number Partitioning 
problem. This algorithm is viewed as a “quantum local search” with matrix elements of 
the Green function G r (r = 1, . . . , n - 1) giving the quantum amplitudes of the transitions 
with different number of spin flips r. Our approach is similar to the analysis of a quantum 
diffusion in a disordered medium with the model of disorder defined by the one- and two- 
point distribution functions P(fi), P r , z (fT)- 

We have shown that the conditional distribution of partition residues P r (ST|ft) in the 
neighborhood of a given string formed by all possible r-bit flips depends on the value of the 
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FIG. 5: Logarithmic plot of C min vs n for randomly generated instances of NPP with 25-bit 
precision numbers. Vertical sets of points indicate results of different trials (~ 100 trials for each 
n, except n=17 with 10 trials). Median values of C mi n are shown with rectangles. Linear fit to the 
logarithmic plot of median values for n between 11 and 17 is shown by the line and gives in Cmin « 
0.55n (Cmin ~ 2 0 8n ) . Very close result is obtained for the linear fit if all data points are used 
instead of the median values. Insert: plot of C(T ) vs T for n=15, precision b=25 bits, do=22. 
Point 1 indicated with the arrow refers to the minimum value of complexity at T = T m j n = 22.67 
where the total population of a ground level po(?min) = 0.15. Point 2 refers to the value of T where 
Po(T) = 0.7. 

partition residue for that string but not on the string itself. This is a specific property of 
the random Number Partitioning problem. 

We used the above property to describe a quantum diffusion in the energy space 
(Eq. (43))- This reduction in the dimensionality leads to the formation of the eigenstate 
which is extended in the energy space. Near the avoided-crossing the adiabatic ground state 
changes from extended to mostly localized near the solution to the optimization problem. 
Because the extended and localized state amplitudes are nearly orthogonal to each other the 
repulsion between the corresponding branches of eigenvalues (the minimum gap) is expo- 
nentially small, g min ~ n2 - "/ 2 , and the run time of the algorithm scales exponentially with 
n. Analytical results are in qualitative agreement with numerical simulations of the time- 
dependent Schrodinger equation for small-to-moderate instances of the Number Partitioning 
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problem ( n < 17). 

One can show that the effect of quantum diffusion in reduced-dimensional space that 
leads to the formation of the extended state can also occur in other random NP-complete 
problems [31]. The method developed in this paper will be applied to study the performance 
of continuous-time quantum algorithms for different random combinatorial optimization 
problems. Also the present framework can be applied to the analysis of quantum annealing 
algorithms for combinatorial optimization problems [32, 33]. This is a classical algorithm 
that is conceptually very close to the quantum adiabatic evolution algorithm considered 
above [34]. The former uses the Quantum Monte Carlo method to simulate on classical 
computers a partition function and ground-state energy of a quantum system with slowly 
varying Hamiltonian that merges at the final moment with the problem Hamiltonian of a 
given classical optimization problem. Among other possible applications of our method is 
the analysis of tunneling phenomenon in the low-temperature dynamics of random magnets. 

We note that the specific property of the Number Partitioning problem (that distinguishes 
it from the other NP-complete problems) is a very weak dependence of P r (ft'|ft) on ft for not 
too large values of ft', ft « y/r(n - r) that takes place for all values of r € [1, n - 1]. This 
rapid fall-off of correlations during the local search (both classical and quantum) is a reason 
that the exponential complexity of optimization algorithms for the Number Partitioning 
problem can be seen already for the relatively small values of n < 15 (cf. Fig. 5). 

Finally, our analysis of sub-harmonic resonances in the Fourier transform I(s) of the 
distribution function P(ft) suggests a possible connection between NPP and the integer 
factorization problem. If, for a given set of a/s, there is a number q that satisfies the 
condition (A3) then dividing all numbers a 3 by q we obtain a new instance of NPP with 
numbers k 3 = a 3 /q that will be completely equivalent to the old one. It is important that 
the precision of the numbers k 3 is restricted by b - log 2 q. If the value of q is sufficiently 
large, log 2 q » b - n, then fc/s correspond to a low precision instance of NPP, i.e. to 
the computationally easy phase mentioned in Sec. II. This is exactly the case when sub- 
harmonic resonances become substantial. One can fix the parameter £ = b/n » 1 in a 
high-precision (computationally hard) case and compute, for randomly generated instances 
{a,} an approximate greatest common divider , i.e. a largest number q that satisfies (A3). The 
distribution of these numbers determines a fraction of high-precision instances of NPP (out 
of all possible 2 nb problem instances) that really belong to a low-precision (computationally 
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easy) “phase”. 

Advance knowledge of this information would be of importance if one is using NPP 
for encryption purposes [17], especially because NPP is otherwise a very difficult problem 
for both quantum and classical computers [29]. It is not obvious at this stage what the 
asymptotic form of this distribution will be in the limit of large n (cf. Fig. A). 

We are not aware of any classical algorithm that could verify if such a number q exists 
for a given set of a 3 in a time polynomial in both n and b. However, on a quantum computer 
one can apply a Shor algorithm to test in polynomial time if strong sub-harmomc resonances 
exist. This question is deferred to a future study. 
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APPENDIX A: SUB-HARMONIC RESONANCES 


We note that function I(s) in (5) can also have additional sharp resonances in the range 
0 < |s| < 2 b . To understand their origin we consider first a particular case when rational 
6-bit numbers cti, a 2 , ■ ■ ■ , o, n all have a number q > 2~ b as a ‘'common divisor’, i.e., there 
exist integers k\, k 2 , ■ ■ ■ ,k n such that 


01 
k 1 


a 2 
k 2 


O'Tl 

k n 


q- 


(Al) 


In this case additional resonances of I(s) occur at the multiples of n jq. Assume now that 
q is no longer an exact divisor of numbers a,- but all the residues of the divisions <Xj/q 
are sufficiently small. Then contributions from the additional resonances at s a rmr/q 
(m = 1, 2, . . .) to the integral in (5) can be estimated as follows (for simplicity we give a 
result for the case E < n 1 / 2 ): 

2” 


P( o) 


y / 27rna 2 (0) 


=- 7 ( 9 ) 


E c (=£*) (-ir 

m= 1 \ / 


(A2) 


p = E 

>= 1 l 


a,* 


7 (?) = 





y^TT ncr 2 (0) 



Here [x] and { 2 :} denote integer and fractional parts of a number x, respectively. If the total 
“dephasing” factor e -7 ^ ~ 1, then contribution (A2) cannot be neglected in the steepest- 
descent analysis of (5) (in general, on should keep contributions from all divisors q with 
small dephasing factors e -7 ^). 


We note that the window function £ 


( m irq \ 

29 ) 


1 for q 2 ” and it decays to zero 


at smaller values of q. We studied numerically the greatest root q max of the the algebraic 
equation 

7 (q) = 7 c ( A3 ) 


for a fixed value of j c < 1. For the sets of random 6-bit numbers dj the dependence of 
the mean value of q min on the problem size n < b is shown in Fig. 6. For n <C 6 we have 
exponential decrease of g max with n and for larger values of n < 6 the value of q m in steeply 
drops to 1. According to the discussion above, in order to neglect the saddle-points with 
s > 0 in (5) (additional resonances) the value of ^ ra i n should satisfy the following condition 
in the asymptotic limit 6 -» 00 : 


<?max < max [2 n , 2 fr ] , 1 < n < 6, 


(A4) 
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with 7 C fixed at some small constant value. Because the precision b that we used in the 
simulations was not very high (limited by machine precision) it is not possible to obtain the 
asymptotic form of the dependence of 5 m j n on n in the range given in (A4). Neither we can 
describe the shape of the plot in Fig. 6 analytically in that range. However, it appears from 
the figure that the condition (A4) is satisfied for sufficiently large n. 

2 25 
2 2 ° 

qmin 2 ' 5 

2 10 

2 s 

2 ° 

FIG. 6: Log-Log plots of the mean value of the largest root of Eq. (A3) <7 m j n vs n. Three sets 
of data points are plotted. Each set of points represents averaging over 25 randomly generated 
instances of NPP. Precision of the random numbers aj is 30 bits and the value of j c = 0.5. Dashed 
line corresponds to the plot of const x 2~ n vs n. Insert: Variance of the log 2 q vs n based on 25 
sample points for each n. Distribution of q m j n values become very broad when the mean drops to 

9min ~ 1- 



APPENDIX B: PROPERTIES OF THE CONDITIONAL DISTRIBUTION OF 
SIGNED RESIDUES IN NPP 


We perform the summation over the spin configurations in Eq. (17) with Eq. (15) taken 
into account. Similar to the derivation of Eq. (5) we use integral representation for delta 
function and obtain 


Prfatf) 





(Bl) 


Uj(s, s') = TTc°s(qj(£ - s')) x JJcos(a,l 
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Here the sum is over all possible subsets J = {juh, ■ ■ ■ ,h ) of length r obtained from the 
set of integers j = 1 , 2 , .... n. Window function <(*) is defined in (5). After the change of 

variables 

x' = s + s', x-s-s\ ( B3 ) 

we obtain from (Bl) that Uj(s, s') factorizes into a product of two terms 

Uj(s,s') = Vj(x)Vj(x') 

^ n «*(% *). Vs(x') = exp ( ) n cos(a, *')■ 


Vj(x) = exp 


je J 


(B4) 


In what follows we will analyze several limiting cases. 
r, n — r f : 

In this case both functions Vj(x) and Vj(x') are very steep and similar to the analysis 
in Sec. II A integrals in (Bl) can be evaluated by the steepest descent method. With the 
appropriate choice of the coarse-graining windows Aft, Aft' in (Bl) (see below) contribution 
to the integrals comes from the vicinity of the point (x = 0,x' = 0). Near this point we use 


PJ cos(aj x) « exp 
je J V 


P[cos(a J , x) « exp I - 
it J 


(n - r)(x' ctj ) 2 
2 


(B5) 


where 

M 2 = ;£4 <^) 2 = — £4 

j€J 

Since each sum here contains a large number of terms we obtain for l.i.d. random numbers 
a u ...,a n (cf. (6)) 


(<jj) 2 « <j 2 (0) + ° (J) > M 2 ~ a2 ^ + ° (^“T 


(B7) 


where a 2 - (a 2 ) is given in (6). Using Eqs. (B4)-(B7) and replacing the window functions 
in (Bl) by unity, we compute the Gaussian integrals in (Bl) and obtain 

v\2 tc\ i rv\2' 


p r (ft,ft') = 


i 


47r<j 2 (0) v /r(n - r) 


exp 


1 /( ft _- ft ') 2 | (ft + ft ') ; 


8<J 2 (0) 


n — r 


(B8) 


The size of the coarse-graining windows in (Bl) is chosen to satisfy the conditions 

2~" j < Aft Aft' < vM n ~ r ) 
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From Eq. (B8) and Eq. (6) one can directly obtain the conditional distribution function 


P T {£l\Q.') 


p r {tf |ft) 


1 

= exp 

y'2 / xna 2 (q) 


(ft'-gft ) 2 

2na 2 (q) 


(B9) 


r = 1; r = n — 1: 

For r = 1 function Vj(x') contains a product of n — 1 terms and is very steep. The 
corresponding integral over x 1 in (Bl) should be taken by the steepest descent method. 
However Vj(x) simply oscillates at frequencies (ft - ft')/ 2 ± % and the integral over x in 
(Bl) should be evaluated using the corresponding oscillating factors. In the opposite case 
r = n - 1, function Vj(x) is very steep and the integral over x in (Bl) should be taken 
by steepest descent. But the integral over x’ there should be evaluated using Vj(x ) that 
oscillates at the frequencies, (ft + ft')/2 ± a r Finally, one can obtain using i.i.d. numbers 
a/s in [0, 1] interval : 

P T (ft'|ft) = i[0(ft T^ + 2) -0(ftTft'-2)] + 0 , (r = l, n-1). (BIO) 


The minus (plus) sign in (BIO) corresponds to r = 1 (r = n - 1). Similarly one can obtain 
the result for any fixed value of r or n — r (that does not scale with n). For |ft|, |ft | < 1 
(BIO) is reduced to (19). 


Numerical simulations of conditional distribution P Tt z (ft ) 


We compute the following integrated quantity: 

i r Q ' 

Q = 2 J drjP r p{ji)i 


(B 1 1) 


for different values of r, ft' and different strings z with E 2 < 1. Numerical results are 
compared in the insert to Fig. 2 with theoretical result below obtained using P r (ft |ft) from 
Eq. (B9) 

(B12) 


i r a> ( ft' \ 

2 1 = 


lo \a(q)\/2nj 

Theoretical and numerical curves nearly coincide with each other. To accurately compare 
the normalization factor in (B9) (see also (20)) we compare the theoretical results with 
numerical values of P r , z (0) for different r and strings z corresponding to E z « 1. The 
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results are plotted in Fig. 2. 
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